#==================================================================================
#   Figures 5a,b,d; 6a,c,d: Calculate and plot mean monthly precip, 
#     community-level SPI index, terms of trade, coffee prices; render other plots
#   Run Stata file paper_data_analysis.do first for input files
#
#   This version: William A. Sundstrom, 11/18/2020  
#==================================================================================

  # clear environment, keep functions
    rm(list = setdiff(ls(), lsf.str())) 

    library(tidyverse)
    library(ggplot2)
    library(SPEI)
    library(haven)
    library(grid)
    library(gridExtra)
    library(ggpubr)

  # Set working directory to root 
    root <- "/Volumes/GoogleDrive/My Drive/Nicaragua Food Security Co-ops Water and Revolution/04 Projects/Paper 3 Nica Livelihoods and Drought 2014-17 /analysis"
    setwd(root)

### Figure 5a
  # read data
    data_5a <- read_dta("output/F5a_data.dta") %>%
      gather(key = type, value = percent, p_short_food:p_short_water) %>% 
      mutate(month = factor(month, levels = month[order(unique(mo))]))
  # plot
    ggplot(data_5a) +
      geom_col(aes(x = month, y = percent, fill = type), position = position_dodge()) +
      ylab("Percent of Households") +
      xlab("") +
      theme_classic() +
      theme(legend.position = "bottom",
            legend.title = element_blank(),
            legend.text = element_text(size = 14),
            axis.title.x = element_blank(),
            axis.title.y = element_text(size = 14),
            axis.text.x = element_text(size = 14),
            axis.text.y = element_text(size = 14),
            panel.grid.major.y = element_line(colour = "grey")) +
      scale_y_continuous(limits = c(0, 1), labels = scales::percent) + 
      scale_fill_manual(values = c("#E69F00", "#4169E1"),
                        labels = c("Food  ", "Water")) +
      geom_line(aes(x = mo, y = p_short_any), size=1)
    ggsave("output/fig_5a.png", width = 8, height = 6, units = "in")
    
### Figure 5d
  # read data
    data_5d <- read_dta("output/F5d_data.dta") %>%
      gather(key = Type, value = Price, maize:beans) %>% 
      mutate(month = factor(month, levels = month[order(unique(mm))]))
  # plot 
    ggplot(data_5d) +
      geom_col(aes(x = month, y = Price, fill = Type), position = position_dodge(), width=0.75) +
      ylab("% Deviation from Monthly Average") +
      xlab("Month") +
      theme_classic() +
      theme(legend.position = "bottom",
            legend.title = element_blank(),
            legend.text = element_text(size = 14),
            axis.title.x = element_blank(),
            axis.text.x = element_text(size = 14),
            axis.text.y = element_text(size = 14),
            panel.grid.major.y = element_line(colour = "grey")) +
      scale_fill_manual(values = c("#8c0738", "#E69F00"),
                        labels = c("Beans  ", "Maize"))
    ggsave("output/fig_5d.png", width = 8, height = 6, units = "in")
    
### Precipitation figures
    
### Read community precip data 
    
    ppt <- read_csv("data/surveycomm2017_monppt.csv") %>% 
      select(-row) %>% 
      rename(ppt = value) %>% 
      mutate(idno = loc_id + 100,
             months_run = (year-1981)*12 + month) %>% 
      arrange(idno,year,month)
    
### Fig 5b seasonality of precip
    
    ppt_mo <- ppt %>% 
      group_by(month) %>% 
      summarize(ppt = mean(ppt)) %>% 
      mutate(month = factor(month, labels = c("Jan","Feb","Mar","Apr","May","Jun",
                                              "Jul","Aug","Sep","Oct","Nov","Dec"))) %>% 
      ggplot(aes(x=month, y=ppt)) +
      geom_line(aes(group = 1), colour="blue", size=1) +
      ylab("Mean precipitation (mm)") +
      theme_classic() +
      theme(legend.position = "bottom",
            legend.title = element_blank(),
            legend.text = element_text(size = 14),
            axis.title.x = element_blank(),
            axis.title.y = element_text(size = 14),
            axis.text.x = element_text(size = 14),
            axis.text.y = element_text(size = 14),
            panel.grid.major.y = element_line(colour = "gray")) 
    ggsave("output/fig_5b.png", width = 8, height = 6, units = "in")    
    
### Fig 6a SPI-3 plot 

  # Function to calculate and retrieve SPI by community using SPEI package
  # i= idno (community index), j = SPI time scale: 1, 3, or 12 months
    spi_get <- function(i,j) { 
      df <- filter(ppt, idno==i) %>% 
        select(ppt) 
      spi_out <- spi(df[,'ppt'], j)
      spi_fit <- data.frame(Y=as.matrix(spi_out$fitted), date=time(spi_out$fitted)) %>% 
        mutate(year = floor(date + 1980),
               month = round((date - floor(date))*12) + 1,
               months_run = round((date-1)*12) + 1) %>% 
        rename(spi = ppt) %>% 
        select(-date)
      return(spi_fit)
    }
    
  # test the function
    spi_get(101,3)

  # loop through and get the SPI for each community at different time frames
    spi_dat2 <- data.frame()
    spi_j <- c(1,3,12)
    indnos <- c(101:152)
    for (jj in spi_j) { 
    for (ii in indnos) { 
     spi_dat1 <- spi_get(ii,jj) 
     spi_dat1 <- spi_dat1 %>% 
       mutate(idno = as.numeric(ii),
              spi_mo = as.numeric(jj))
     spi_dat2 <- bind_rows(spi_dat2, spi_dat1)
     } }
   
   spi_comm <- spi_dat2
   
### Plot SPI for 2010-2017 with shaded May-Oct. (rainy months)
   
    spi_plot <- spi_comm %>% 
     filter(year>=2010 & year<=2017) %>% 
     group_by(months_run, spi_mo) %>% 
     mutate(med_spi = median(spi)) %>% 
     ungroup()
    May <- spi_plot %>% 
     filter(idno==101 & spi_mo==3 & month==5) %>% 
     mutate(May = months_run) %>% 
     select(year,May) 
    Oct <- spi_plot %>% 
     filter(idno==101 & spi_mo==3 & month==10) %>% 
     mutate(Oct = months_run) %>% 
     select(year,Oct) 
    rainy <- left_join(May,Oct)
   
  # function for plots: j is SPI horizon = 1,3,12
    spi_plotter <- function(j) { 
     yt <- paste0("SPI-",j)
     plot <- spi_plot %>% filter(spi_mo==j) %>% 
       ggplot() + 
       geom_rect(data=rainy,aes(xmin=May, xmax=Oct, ymin=-Inf, ymax=+Inf), fill='pink', alpha=0.4) +
       geom_line(aes(x=months_run, y=spi, group = as.factor(idno)), color="light blue", size=0.2) + 
       geom_line(aes(x=months_run, y=med_spi), color="dark blue") + 
       geom_vline(xintercept=c(402,438), color="black", size=1, linetype="dotted") + 
       labs(y = yt) +
       theme_classic() +  
       theme(axis.title.x=element_blank(),
             panel.grid.major.y = element_line(colour = "grey"),
             panel.grid.major.x = element_line(colour = "grey"),
             axis.text.x = element_text(size=10),
             axis.text.y = element_text(size=10)) + 
       scale_x_continuous(breaks = c(349,361,373,385,397,409,421,433,445,457),
                          labels = c("Jan 2010","Jan 2011","Jan 2012","Jan 2013","Jan 2014",
                                     "Jan 2015","Jan 2016","Jan 2017","Jan 2018","Jan 2019")) +
       scale_y_continuous(breaks = c(-3,-2,-1,0,1, 2, 3, 4), limits = c(-3,4)) 
     return(plot)
    } 
   
    spi_plotter(3)
    ggsave("output/fig_6a.jpg", width = 7.5, height = 3.0, units = "in")
   
### Fig 6c Terms of trade: Corn and beans vs coffee 
   
    fig6c <- as.data.frame(read_dta("output/F6cd_data.dta")) %>% 
      filter(yy>=2010 & yy<2018) 

    f6c <- fig6c %>% 
         ggplot() + 
         geom_line(aes(x=mt, y=tt_beansred_coffee_jf), color="maroon", size=.75) + 
         geom_line(aes(x=mt, y=tt_maize_coffee_jf), color="#E69F00", size=.75) + 
         geom_vline(xintercept=c(654, 690), color="black", size=1, linetype="dotted") + 
         theme_classic() +  
         theme(axis.title.x=element_blank(), axis.title.y=element_blank(),
               panel.grid.major.y = element_line(colour = "grey"),
               panel.grid.major.x = element_line(colour = "grey"),
               axis.text.x = element_text(size=10),
               axis.text.y = element_text(size=10)) + 
         scale_x_continuous(breaks = c(600,612,624,636,648,660,672,684,696),
                            labels = c("Jan 2010","Jan 2011","Jan 2012","Jan 2013",
                                       "Jan 2014","Jan 2015","Jan 2016","Jan 2017","Jan 2018")) + 
         annotate("text", x = 618, y = 260, label = "Red beans") + 
         annotate("text", x = 680, y = 220, label = "Maize")
    f6c
    ggsave("output/fig_6c.jpg", width = 7.5, height = 3.0, units = "in")
    
### Fig 6d Coffee prices 
    
    fig6d <- as.data.frame(read_dta("output/F6cd_data.dta")) %>% 
      filter(yy>=1990 & yy<2018) 
    
    f6d <- fig6d %>% 
      ggplot() + 
      geom_line(aes(x=mt, y=OtherMilds), color="dark green", size=.75) + 
      theme_classic() +  
      theme(axis.title.x=element_blank(), 
            panel.grid.major.y = element_line(colour = "grey"),
            panel.grid.major.x = element_line(colour = "grey"),
            axis.text.x = element_text(size=10),
            axis.text.y = element_text(size=10)) + 
      scale_x_continuous(breaks = c(360,420,480,540,600,660),
                         labels = c("Jan 1990","Jan 1995","Jan 2000","Jan 2005",
                                    "Jan 2010","Jan 2015")) +
      scale_y_continuous(breaks = c(0,50,100,150,200,250,300), limits = c(0,310)) + 
      ylab("U.S. cents per pound")
    f6d
    ggsave("output/fig_6d.jpg", width = 7.5, height = 3.0, units = "in")
    
